## ---------------------------------------------------------------------- #
## Figure 6
## Christopher J. Fariss
## "Yes, Human Rights Practices Are Improving Over Time"
## American Political Science Review
## ---------------------------------------------------------------------- #
rm(list = ls())

library(rstan)

load("Full_Models_PPP/M_All.RData")
ppp.all.out <- extract(mod)
rm(mod)

load("Full_Models_PPP/M_Fixed.RData")
ppp.fix.out <- extract(mod)
rm(mod)

load("Full_Models_PPP/M_Standards.RData")
ppp.stand.out <- extract(mod)
rm(mod)


load("Full_Models/M_All.RData")
theta.all.out <- extract(mod)
rm(mod)

load("Full_Models/M_Fixed.RData")
theta.fix.out <- extract(mod)
rm(mod)

load("Full_Models/M_Standards.RData")
theta.stand.out <- extract(mod)
rm(mod)


NAMES <- c("PTS:\nState",
"PTS:\nAmnesty",
"PTS:\nHuman Rights Watch",
"Hathaway:\nTorture",
"ITT:\nTorture",
"CIRI:\nTorture",
"CIRI:\nExtrajudicial Killing",
"CIRI:\nDisappearance",
"CIRI:\nPolitical Imprisonment",
"Harff:\nGenocide",
"Ufelder and Valentio:\nMass Killing",
"Harff and Gurr:\nMassive Repression",
"Rummel:\nGenocide and Politicide",
"UCDP:\nOne-Sided Killing",
"WHPSI:\nNegative Sanctions",
"WHPSI:\nExecutions"
)

data.NAMES <- list("State",
"Amnesty",
"HRW",
"hathaway",
"ITT",
"TORT",
"KILL",
"DISAP",
"POLPRIS",
"genocide",
"mass_killing",
"massive_repression",
"rummel",
"killing_present",
"negative_sanctions",
"executions"
)



cor.theta.stand <- cor.theta.fix <- cor.theta.all <- list()

par(mfrow=c(4,4), mar=c(2,2,1,1))

for(k in 1:16){
    for(j in 1:3){
        if(j==1) data <- read.csv("M_2_Full_Data.csv")
        if(j==2) data <- read.csv("M_1_Full_Data.csv")
        if(j==3) data <- read.csv("M_All_Vary_Full_Data.csv")
        
        SIM <- 1000
        mat <- matrix(rnorm(matrix(1,nrow=nrow(data),ncol=SIM), mean=data$theta_mean, sd=data$theta_sd), nrow=nrow(data),ncol=SIM)
        
        
        YEARS <- unique(data$YEAR)
        
        
        temp <- data[,paste(c("COW", "YEAR", "id", data.NAMES[[k]]))]
        
        
        observed.year.means <- lapply(1:length(YEARS), function(i){
            mean(temp[temp$YEAR==YEARS[i],data.NAMES[[k]]], na.rm=T)
        })
        
        
        theta.year.means <- lapply(1:length(YEARS), function(i){
            apply(mat[temp$YEAR==YEARS[i],], 2, mean)
        })
        
        
        
        Mu <- unlist(lapply(theta.year.means, mean))
        
        
        mat2 <- do.call("rbind", theta.year.means)
        
        test <- lapply(1:SIM, function(i){
            cor(unlist(observed.year.means), mat2[,i], use="pairwise", method="spearman")
        })
        
        if(j==1)cor.theta.stand[[k]] <- unlist(test)
        if(j==2)cor.theta.fix[[k]] <- unlist(test)
        if(j==3)cor.theta.all[[k]] <- unlist(test)
        
        
    }
    #boxplot(cor.theta.all[[k]], cor.theta.fix[[k]], cor.theta.stand[[k]])
    
}


#par(mfrow=c(4,4), mar=c(1,3,3,.5));for(k in 1:16){boxplot(abs(cor.theta.stand[[k]]) - abs(cor.theta.fix[[k]]), las=2, ylim=c(-.8,.8)); abline(h=0); mtext(side=3, NAMES[k])}


#par(mfrow=c(4,4), mar=c(1,3,3,.5));for(k in 1:16){boxplot(abs(cor.theta.stand[[k]]) - abs(cor.theta.all[[k]]), las=2, ylim=c(-.8,.8)); abline(h=0); mtext(side=3, NAMES[k])}


par(mfrow=c(4,4), mar=c(.5,4,4,.5))
for(k in 1:16){
    boxplot(abs(cor.theta.stand[[k]]) - abs(cor.theta.all[[k]]), abs(cor.theta.stand[[k]]) - abs(cor.theta.fix[[k]]), las=2, ylim=c(-1,1), xaxt="n", lwd=.5)
    abline(h=0, col=grey(.5), lwd=.5)
    #abline(v=1.5, col=grey(.5), lwd=.5)
    mtext(side=3, NAMES[k], line=.5)
    
    mtext(side=3, expression(italic("Standard")), line=-1, at=1.5, cex=.7, font=1)
    mtext(side=1, expression(italic("All Varying")), line=-1, at=.95, cex=.7, font=1)
    mtext(side=1, expression(italic("Constant")), line=-1, at=2.05, cex=.7, font=1)
    
}

